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We report on a numerical simulation of the classical evolution of the plane-wave matrix model with 
semiclassical initial conditions. Some of these initial conditions thermalize and are dual to a black 
hole forming from the collision of D-branes in the plane-wave geometry. In particular, we consider 
a large fuzzy sphere (a D2-brane) plus a single eigenvalue (a DO-particle) going exactly through 
the center of the fuzzy sphere and aimed to intersect it. Including quantum fluctuations of the 
off-diagonal modes in the initial conditions, with sufficient kinetic energy the configuration collapses 
to a small size. We also find evidence for fast thermalization: rapidly decaying autocorrelation 
functions at late times with respect to the natural time scale of the system. 



- Introduction and conclusion. From the beginning of 
the gauge/gravity correspondence [T] it was understood 
that large black holes with anti de Sitter asymptotics 
should be related to thermal states in the dual field the- 
ory [5]. So, formation of a black hole from non-thermal 
initial conditions should be dual to thermalization of a 
specific initial state in the dual dynamics. This idea has 
led to the suggestion that the rapid thermalization ob- 
served in heavy ion collisions |3] should be described by 
a dual black hole formation event 

Nascent black holes settle very quickly to the no-hair 
solutions, as shown by numerical simulations (see, e.g., 
[5]). We expect that the dual theory will behave simi- 
larly. It has also been conjectured that black holes are 
"fast scramblers" [6], i.e., they distribute information (or 
wash it away) faster than any other physical system, log- 
arithmically in the number of degrees of freedom. This 
refers to how fast small quantum fluctuations away from 
equilibrium settle back to equilibrium in a black hole. 

The purpose of this paper is to explore such thermal- 
ization processes from the point of view of the field theory 
and to formulate a program where the fast scrambler con- 
jecture can eventually be tested by numerical methods. 
We do this by focusing on a system with finitely many 
degrees of freedom where a collision of two gravitons or 
D-branes at high energy can in principle be studied: the 
plane- wave (or BMN) matrix model [7]. This paper de- 
scribes the first simulations of this system we have per- 
formed and the evidence we have acquired for fast ther- 
malization. The simulations solve the classical equations 
of motion of the model. We show that the time averages 
of quantities in the system after thermalization match 
the Gibbs distribution for some degrees of freedom that 
appear quadratically in the Hamiltonian. The temper- 
atures measured by various of these degrees of freedom 
are the same. We also show that various gauge invariant 
quantities have autocorrelation functions rapidly decay- 
ing in time with respect to the natural time scale in the 
system (defined by the data rather than machine time) . 

Ideally we would study this problem in the matrix 



model of [8], which is dual to M-theory on fiat space 
as the discrete light-cone formulation of the theory. In 
that system one has asymptotic states that can be scat- 
tered and could lead to a black hole formation event with 
an S-matrix interpretation. However, the initial condi- 
tions for that setup are not understood: the gravitons 
are bound states at threshold that necessitate solving the 
many body quantum dynamics in detail. Monte-Carlo 
simulations of the matrix dual black holes of this model 
have shown that the ensemble is unstable due to the flat 
directions of the matrix model potential . The instabil- 
ity can be regulated with mass terms that are naturally 
present in the BMN model. Euclidean computations of 
the BMN model at equilibrium have been performed in 
jlOj . This paper deals with the time-dependent classical 
regime in the theory and the dynamics of thermalization. 

The BMN matrix model, which represents the plane- 
wave geometry with maximal supersymmetry, solves the 
problem of describing gravitons by having a different clas- 
sical solution for each graviton. Each of these gravitons 
is represented by a fuzzy sphere. We lose the ability 
to perform collisions with asymptotic states that scatter 
from each other. However, as shown in |11[ . there are 
initial conditions where the fuzzy spheres are displaced 
with respect to each other and we can set up periodic 
brane collisions (crossings) instead. The period for these 
setups, r, is independent of the details of the graviton 
branes and their velocities. This lets us measure time 
with a clock that has a well-defined physical interpreta- 
tion. These crossings have classical instabilities associ- 
ated with them: some degrees of freedom are tachyonic 
during a brief time of the r-period and grow exponen- 
tially in the number of crossings until the system back- 
reacts. The growth of fiuctuations has been understood 
with a linearized analysis [llj and in the present letter 
we extend that analysis to the rest of the thermalization 
process where the dynamics is very non-linear. Classi- 
cally these fluctuations can be set to zero, but they will 
be present in the full dynamics because of quantum me- 
chanics. 
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We add noise to represent quantum fluctuations in the 
initial conditions. Our configurations seed the modes 
that would become tachyonic during some portion of the 
T-period . We initialize each such mode with a gaus- 
sian probability function with a width determined by h 
and the adiabatic frequency of the mode. Once these fluc- 
tuations grow sufficiently they take over and scramble the 
system. The time scale for the initial exponential growth 
is logarithmic in the size of the initial fluctuations, which 
are proportional to \/h. We consider the system to ther- 
malize quickly if, subsequent to this period, the decay of 
fluctuations back to equilibrium is fast. The analysis we 
do is valid strictly only when h is small. The solutions 
have large energy of order N'^, where N is the size of the 
matrices. The energy does not scale with h. If these sys- 
tems thermalize, their temperature is large in quantum 
units. This is exactly the regime where classical physics 
is valid, for we only have finitely many degrees of free- 
dom and all of them will have large quantum numbers 
and can be treated classically (this sets the limit of how 
big we can make h in practice). This means in particular 
that we can ignore the fermions, for they only affect the 
low temperature regime. 

In the rest of this letter we discuss the numerical im- 
plementation of the BMN matrix model and the evidence 
for fast thermalization of the initial conditions we chose. 

- Numerical implementation. The bosonic degrees of 
freedom of the BMN matrix model are the hermitian ma- 
trices X'^^^'^'^ and y^=i'- -.6 Qj^^ their canonical conju- 
gates Pi and Qa- The bosonic part of the Hamiltonian 
is 

H =^tr(^Pl + Ql + {X' + ie'i^X^X^f 

+i(y")2 - [x'^vY - ^yy'^y^?)- (1) 

We have rescaled the variables so that the classical equa- 
tions of motion are independent of h and all the quantum 
mechanics is hidden in the initial conditions. We have 
also normalized the mass of X to one, i.e., we measure 
time by the oscillation period of one of the X modes. 

Because of the U(A^) gauge symmetry we must enforce 
the Gauss' law constraint: 

C =[X\P,] + [Y'^,Qa]^0. 

To solve the equations of motion we use a leapfrog algo- 
rithm and we record the absolute value of the constraint 
tr(C^) as a check for the code. We find that the con- 
straint is well satisfied for the runs we perform, so we do 
not need to implement constraint damping. 

The main sources of difficulty are the initial conditions. 
For this paper, we have used the following initial classical 



configuration: 
X 

^0 



-\0 oj'^ ~\5x\ j' 

l,2^0^gl,...,6^ 



X' 



Ll 5x2 
5x1 



V 



P' 



The dimension of the matrices above is set to iV = n + 1. 
The U are SU(2) angular momentum matrices in the 
ri-dimensional representation. This is a fuzzy sphere of 
size n. The system has an additional eigenvalue that is 
initially at the origin with velocity v in the positive X^ 
direction. These are the initial conditions discussed in 
[TTj with the addition of fluctuation seeds 5x, 5y. The 
(5x, 5y are generated randomly using a complex gaussian 
distribution with a width proportional to ^Jhjn. We in- 
terpret these as quantum fluctuations of the off-diagonal 
degrees of freedom. 

Recall that the ground state of an oscillator with 
Hamiltonian IR = p^ + us^x^ has a gaussian wave func- 
tion with squared width (x^) = h/2ui. In our case, 
because of the initial conditions, all of the off-diagonal 
modes between the lone eigenvalue and the fuzzy sphere 
have approximately the same frequency of oscillation, 
proportional to n [12] . 

The 5x^ are determined by forcing the matrices to be 
hermitian. All the off-diagonal 5y components are gener- 
ated by the same gaussian distribution, whereas the diag- 
onal ones are set to zero since they are subleading in N . 
This is a very rough approximation for the Y modes con- 
necting the fuzzy sphere to itself, lumping them together 
as if they all had the same mass. We do not add fluctu- 
ations in the modes connecting the fuzzy sphere to itself 
in the X variables as the unstable modes grow so quickly 
that such fluctuations are not required. If the system 
thermalizes the fine details of the initial conditions get 
washed out at later times, so we only need them to be 
qualitatively correct. Notice that our initial conditions 
are built to exactly satisfy C — while preserving the 
typical size of quantum fluctuations for the space vari- 
ables. This is why we have no fluctuations of the P, Q 
variables nor of the off-diagonal modes of X^ connecting 
the lone eigenvalue and the fuzzy sphere. 

The discretized matrix equations of motion read 

dV 
' dX 

and similarly for the Y modes. Here V is the poten- 
tial obtained from eq. ([T]). The parameter St and the 
total number of iterations of the leapfrog algorithm can 
be varied in the numerical code. We record the matrix 
configurations every few steps. 

Due to the accumulation of numerical rounding errors, 
every few steps we need to force the matrices to be her- 
mitian. We do this right before recording configurations. 
We also vary the random seed to generate an ensemble 
with gaussian distributions and check that the results are 



t-rPt+itSt, Pt+it=Pt^ 



6t. 
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robust against these variations. A more detailed descrip- 
tion of the code and a more comprehensive presentation 
of the numerical results and their interpretation will be 
presented elsewhere (also see supplement). 

- Results. We now describe some useful ways to visu- 
alize the information contained in the X, Y modes and 
their time derivatives. To describe the thermodynamics, 
we need to coarse grain the degrees of freedom, which 
must be gauge invariant combinations of X and Y. The 
simplest such combinations are traces of matrix products. 
We can use these to compare different values of N and to 
study the large N thermodynamic limit. Note that the 
traces of the matrices X* and Y'^ are decoupled: the non- 
linear parts of the equations of motion are commutators, 
so the traces of the matrices evolve independently from 
the rest of the system. The trace of X* oscillates with 
angular frequency w = 1, while the trace of Y"' oscillates 
with Lu = 1/2. Because of our initial conditions, only the 
trace of the X° mode is excited and it serves as our clock. 

The other invariant way to work with matrices is in 
terms of their eigenvalues. These can tell us about the 
dual D-brane geometry. When the matrices are approxi- 
mately commuting there is a clear geometric interpreta- 
tion: the eigenvalues are positions of D-branes. When 
the matrices do not commute they still serve to roughly 
describe the distribution of D-branes inside the fuzzy ob- 
ject. In Fig. [l] we plot the eigenvalues of X^(t). Initially 



dard deviations of the eigenvalues of the matrices, see 
Fig. [2] As shown in the figure, the fuzzy sphere collapses 




FIG. 1. Eigenvalue evolution for X". Here and in the fol- 
lowing figures we have set n = 10, v — 20, and h = 0.001. 
The time axis is in discrete time units between recordings of 
configurations. 

all of the motion is in the lone eigenvalue and the other 
eigenvalues are evenly spaced - a property of the fuzzy 
sphere. As time goes by, the eigenvalues collapse to a 
much smaller vertical extent and oscillate collectively. 
The eigenvalues become very unevenly spaced and upon 
zooming in appear to repel each other, showing a typical 
behavior of random matrices. Qualitatively, they behave 
like a Dyson gas [13] , but a detailed comparison is beyond 
the scope of the present paper. 

It is also interesting to study the size of the system in 
different directions. We do this by evaluating the stan- 
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FIG. 2. Standard deviation of eigenvalues for the matrices 
X°'^ and We use the trace of rescaled, to keep 

track of time (black curve at the bottom). Other values of 
n, V, and h are qualitatively similar, see supplement. 

in size substantially. After the sphere has largely col- 
lapsed, the Y modes grow from zero and converge to a 
value that is very close to the late-time value for the X 
modes. Their growth is controlled by the random time 
variation of their effective mass after collapse. It is be- 
cause of this that we needed to include fluctuations for 
most of the Y modes. The subset of Y modes connecting 
the fuzzy sphere and the eigenvalue do not grow enough 
in the initial phase and the time for fluctuations of those 
modes to converge is substantially longer. The size has 
only small fluctuations after convergence and the system 
seems to stabilize rapidly. The figure suggests that the 
object is becoming nearly spherical, a property shared 
by black holes without angular momentum. However, 
the corresponding dual black holes should have some de- 
formation since they are not in asymptotically flat space. 

To test for thermalization, we compare time averaged 
distributions over successive configurations to those of 
the Gibbs ensemble for the classical system at some tem- 
perature T. Using the Gibbs measure dP dQ exp{—H/T) 
we see that the momentum variables factorize into gaus- 
sian integrals. Thus the momenta are determined by 
the gaussian ensemble for hermitian matrices. It is well 
known that the distribution of eigenvalues (sufficiently 
coarse grained) should be a semicircle. We test this for 
the P° and matrices starting way after the system 
looks thermalized {e.g., after t ~ 600 in Fig. [2|. We 
wait until t = 5000 to measure thermal properties just to 
make sure. This is shown in Fig. [3] The semicircle model 
matches the data well for both the X and Y momenta, 
which have the same distribution. This suggests that the 
system has thermalized, as the temperature measured 
from the AT's is the same as that measured from the F's. 

Now that we have numerical evidence for thermaliza- 
tion we can study near-equilibrium configurations and 



FIG. 3. The eigenvalues of P° and are binned from 
t = 5000 to t = 20000 every ten steps, using the same time 
units of Fig. [2] The semicircle distribution is integrated over 
the binning intervals and normalized to the total count of 
eigenvalues. The width is matched using the standard devia- 
tion of the eigenvalue distribution. 



fluctuation decay rates. This information can be ob- 
tained from the autocorrelation function {0{t)0^{t + a)), 
where 0(t) is some classical gauge invariant observable. 
This is an application of the fluctuation-dissipation the- 
orem. This is averaged over t well after thermalization. 
The simplest observables we can consider are of the form 
tr{X^ +iX'^)'" . Similar traces are identified with graviton 
modes in A/" = 4 super Yang-Mills theory [T] , where they 
are interpreted as having angular momentum L along the 
dual . Here, L denotes angular momentum in the 12 
plane of the X variables. Higher L values correspond 
to higher spherical harmonics in the dual geometry. For 
i = 1, as we have seen, the mode is decoupled, so the 
simplest non-trivial case will be for L — 2, shown in Fig. [4] 
(top) . We can see that the autocorrelation for L = 2 dies 
off quickly, with respect to the natural external clock. We 
also note that it oscillates in an interesting pattern, indi- 
cating that there are internal oscillation times of the vari- 
ables associated with the thermalized system. Relative 
to these internal oscillations the autocorrelation function 
decays quickly (by half within 2 oscillations), so the as- 
sociated vibration modes have a low quality factor. This 
indicates fast thermalization. In Fig. |4] (bottom) we com- 
pare autocorrelations for higher L. We notice that the 
higher the L, the faster the autocorrelations decay. This 
is expected from black hole physics and the membrane 
paradigm of the horizon: when information approaches 
the membrane, it diffuses along the membrane until it 
becomes uniform. Diffusion happens first at short dis- 
tance scales and then cascades to large scales. So this 
is evidence for an approximate notion of locality in the 
angular directions even in the thermal regime. 
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FIG. 4. Top: normalized autocorrelation of tr(X^ + iX^)^ . 
The width of the band indicates our statistical uncertainty 
from few similar sequences related to each other by rotations. 
We include our clock measure (black curve at the bottom). 
Bottom: we compare different L modes for a shorter time 
period. The configurations are averaged in time starting at 
iteration 5000. Other values of n, v, and h are qualitatively 
similar, see supplement. 
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SUPPLEMENT 

This is a supplement to "Evidence for fast thermaliza- 
tion in the plane-wave matrix model" , as appeared on 
Physical Review Letters. The supplement contains tech- 
nical details on how the data was extracted from the 
simulation, as well as how the data was reduced. We in- 
clude additional plots on studies that show that the plots 
contained in the letter reflect the data broadly. We also 
include some additional discussion on the meaning of dif- 
ferent thermalization times and what are the systematics 
that affects their numerical evaluation. 



Methodology 

- Equations of motion and numerical parameter 
choices. Given the Hamiltonian 

H=^ tr(^if + Ql + (X' + ie'^'^X^X'^f 

+i(y'')2-[x\y«]2-i[r",r'']2), 

the equations of motion for the X and Y fields are dif- 
ferent. They are given by 

X'^P,, 

P, = -X' ~ 3iJ2^'^''[X\X''] 

jk 

k a 

Qa = -\y'' + Y.^[X\Y%X^] + Y.^[Y\Y%Y^] . 

i 3 

Notice that for large values of X, Y and for X, Y some- 
what random the frequency of the modes is parametrized 
roughly by the absolute value of the eigenvalues of X as 
a matrix. Also notice that if we take traces we find that 

tri:^ = trP, , 
tri"^ = trQ, , 
trPi = -tiX' , 

trQ, = -^trr^ 

These are closed sets of linear equations and they are 
harmonic oscillators of frequencies 1 and ^ respectively. 
These give us our clock. The amplitude of irX has a 
periodic motion and this period is independent of the 
details of the initial condition. Checking that this mode 
has this periodic motion serves as a test of the code. 

From this information the time steps between itera- 
tions on our code should be controlled by the maximum 
size of X. For our simulations we have that X is of 
order 10 or 20, so the time steps 5t should be such that 
i5t ^ 27r/20 ~ 0.3 so that a single oscillation of the fastest 



mode is covered by many points. For the particular sim- 
ulation that is quoted in the letter we choose 5t = 0.004 
and we save the configurations every 25 iterations. We 
have done some other configurations with a finer time 
step, but the quality of the results are comparable to the 
ones shown in the letter. 

In this supplement we also consider other values of N 
(in the letter we had considered the 11 x 11 case only). 

- Computing autocorrelation functions. In the letter 
we state that we compute the autocorrelation functions 

fo{a)^{0{t)0\t + a)), 

where 0{t) is some classical gauge invariant observable. 
These are averages over time, so the variable t is averaged 
over. These are good measurements of the dynamics at 
equilibrium. The initial time cutoff is controlled by ask- 
ing when has the system thermalized. In the letter in 
Fig. 2 we have shown that the size of the system reaches 
equilibrium values at some t close to 500 (in iterations 
that are saved). We wait until machine time is 5000 be- 
fore we start averaging. Our simulation runs for 20000 
stored iterations. We get the autocorrelation functions 
from computing the Fourier transform of the power spec- 
trum of the last 15000 saved iterations. These depend on 
the time displacement a and on the observable. 

Due to rotational symmetry of the configurations there 
are various autocorrelation functions that are the same. 
We can use this to get a measurement of statistical error 
bars with a sample of 3 from a single run. We normalize 
these to the standard deviation of fluctuations, so the 
plot depicts the quantities 

/o(«)//o(0) 

and we use the average value /ci(0) to do this. Consid- 
ering we have various data sets, they display some small 
variation in the size of fluctuations. This is expected. 



Additional plots 

- Varying the initial velocity of the single eigenvalue. 
In Fig. [5] of this supplement we plot the standard devi- 
ations of the eigenvalues of the matrices X' and as 
in Fig. 2 of the letter, but for different initial values of 
the velocity v. We go from v — IQ to v = 100 in steps 
of 10. From these plots we see that the confluence time 
at which all the standard deviations of the various vari- 
ables converge to the same value is roughly independent 
of V. This confluence time is approximately given by 15 
- 18 periods of the internal clock (the black curve at the 
bottom of the plots) . 

The rest of the initial conditions are as follows: 

TV = 10 -1-1, /l=10"^ Ai = 0.0008, Ncycies = '2.h0. 
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The variable Ncydes dictates how long we wait before 
writing new configurations to a file. Notice the relatively 
high value for Ncydes ■ We had selected this value because 
we wanted to look at long simulations without generat- 
ing huge dump files, but of course this makes the curves 
look more discontinuous than what we had in the letter 
(especially for X°, the red curve in the plots). 

These also were obtained at a smaller value of h than 
what is reported on the letter. We will have further com- 
ments on this later on in this supplement. These were 
generated in a previous iteration of the code where the 
definition of h is slightly different, thus the tilde on top of 
it. The main reason for the new definition was to make 
it possible to fix the quantity h to be able to compare 
different values of N homogeneously. 

~ Varying the size of the matrices. To show that our 
results are fairly independent of we now include plots 
for = 20 + 1 similar to the ones contained in the letter 
for the TV = 10 -I- 1 case. 

The figures shown in Fig. [6] are calculated with N = 
20 -f 1, w = 60, /i = 10"'', St = 10"''. We record the 
configuration every 20 iterations. The cutoff for ther- 
malization is taken at saved iteration 5000, immediately 
when the size lines coincide. This simulation collected 
25000 iterations from the beginning. The plots here are 
generated in the same way that the same plots used in 
the letter in Figs. 2 and 4 (the top one) we generated. 
They are the plots for the new configuration. The plots 
are qualitatively similar to the ones included in the letter. 
A detailed comparison between different N is something 
we are looking into in order to characterize the system 
better. Notice that the autocorrelation functions have a 
different natural oscillation time relative to the clock. 

In Fig. [7] we also include plots equivalent to the ones 
in Fig. [5] this time keeping u = 40 fixed but changing 
N. Again, we see that the confluence times are quite 
independent of N. 

Thermalization times 

A final technical note that we include here is regard- 
ing the notion of thermalization times. There are three 
notions of thermalization times that we can consider. 

First, we can ask what is the time it takes for a small 
quantum fluctuation (adding a particular quantum bit) 
to scramble so that we can not recover that bit any longer 
as a function of the size of the system. This is the ques- 
tion that is asked in the paper by Sekino and Susskind 
[6], where it is argued that black holes are fast scram- 
blers. This question is for a near equilibrium situation. 
Our code does not provide yet for a way to add additional 
fluctuations after a black hole is formed, so we can not 
address this question directly. Also, we have not explored 
how to change the size of the dual black holes systemat- 
ically: we need to change N keeping something fixed as 



we vary N. It is not clear from the paper [B] what is the 
correct way to implement this idea in our simulations. 

Secondly, there is an intrinsic notion of scrambling near 
equilibrium as provided by the fluctuation dissipation 
theorem. Understanding how the autocorrelations die 
gives us a classical notion of how long it takes for a fluc- 
tuation in a particular channel to die. As described in 
the letter in Fig. 4, and in the supplement in Fig. |6j the 
autocorrelation functions both decay and oscillate. The 
decay of the signal as measured from a = to the next 
peak in the oscillation is about 0.5 — 0.6 and within about 
3 oscillations it has decayed to about less than a third of 
the size. This means that the autocorrelation function 
dies off as quickly as it oscillates and can therefore be 
considered a fast scrambler from this point of view. This 
ratio of times does not seem to depend too strongly on 
the size of the system. 

Third, there is the time to thermalization from a par- 
ticular initial condition. Since the only place where h is 
used in the letter is in the details of the initial condi- 
tions, we can ask what is the h dependence of the time 
to thermalization fixing everything else. When H is zero, 
the configuration is in periodic motion forever and it does 
not thermalize. When h is small, the off-diagonal fluctu- 
ations are tiny and we can solve the equations for these 
modes in the linearized regime. This was done in detail 
in the previous work by some of the authors [TT]. The 
fluctuations are either stable or unstable. The unsta- 
ble modes grow exponentially between oscillations until 
back reaction of the nonlinearities kicks in. The non- 
linearities can be understood as fluctuations reaching a 
particular size. We can call this size of fluctuation the 
saturation size and the time at which the fluctuations get 
to this size the saturation time. Before this, the growth 
is exponential, so the time for this growth depends on 
the size of the initial fluctuation logarithmically. This is 
t cx log{ A gat /Ain). Since Ai^ oc \/^, the time to satura- 
tion grows logarithmically in h. Evidence for this can be 
seen by comparing the second figure in this supplement 
with respect to the corresponding Fig. 2 in the letter. 
Notice also that after the X modes shrink, it takes longer 
for the Y modes to saturate in the plot in this supple- 
ment relative to the one appearing in the letter. This is 
because the Y modes grow slower in the initial stage, so 
their amplitude is small after the first collapse. 

Notice that if we keep X random, but Y infinitesimal 
in the equations of motion, the Y motion is linear in Y . 
So again, if the initial size after the first collapse of X is 
controlled by h, the time to saturation of the Y modes is 
logarithmic in h. 

The plots give evidence for this behavior. Obviously, 
if we increase h we shorten the time to thermalization 
from the initial condition. Care has to be taken on what 
exactly is an equilibrium configuration, so this time de- 
pends on the systematics of how we declare that a par- 
ticular configuration has reached equilibrium. 




FIG. 5. Simulations with different values of the initial velocity v. The time to thermalization is found to be roughly independent 
of V. 
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FIG. 6. Left: simulation of size of configuration in various directions as measured with standard deviation of eigenvalues at time 
t. We include our clock measure (black curve at the bottom). Right: normalized autocorrelation functions of tr[(Xi + iX2)^] 
and comparison to the clock. The width of the band indicates our statistical uncertainty from few similar sequences related to 
each other by rotations. 
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FIG. 7. Simulations with different values of the matrix size N — n + 1. The initial velocity is fixed at u = 40 and the rest of 
the parameters are as in S. The time to thermalization is found to be roughly independent of n (at t ~ 500). 



